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^ • 1 Introduction 

c/2 , 

■ There is substantial evidence that the ground state in many of oxides is 

' inhomogeneousT. In cuprates, for example neutron-scattering experiments 

suggest that phase segregation takes place in the form of stripes or short 
I segments of stripes [H [3] ■ There is some controversy whether this phase segre- 

^ i gation is associated with magnetic interactions. On the other hand it is also 

O I generally accepted that the charge density in cuprates is not homogeneous. 

The idea of charge segregation has a quite long history (see for example [H 
[5j[6j). In charged systems the phase separation is often accompanied by the 
charge segregation. Breaking of the charge neutrality leads to the appearance 
of the electric field and substantial contribution of the electrostatic energy to 
the thermodynamic potential[71 [H]. Recently it was suggested that interplay 
of the short range lattice attraction and the long-range Coulomb repulsion 
between charge carriers could lead to the formation of short metallic [9l [10] 
or insulating [lOl [TT| stripes of polarons. If the attractive potential is isotropic, 
. charged bubbles have a spherical shape. Kugel and Khomskii [12^ suggested 

recently that the anisotropic attraction forces caused by Jahn- Teller centers 
^ . could lead to the phase segregation in the form of stripes. The long-range 

anisotropic attraction forces appear as the solution of the full set of elasticity 
equations (see ref.|l3|). Alternative approach to take into account elasticity 
' potentials was proposed in ref.jilT and is based on the proper consideration 

C , of compatibility constraint caused by absence of a dislocation in the solid. 

Phenomenological aspects of the phase separation was discussed recently in 
the model of Coulomb frustrated phase transition [151 [HI [H] ■ 

Here we consider some aspects of the phase separation associated with 
different types of ordering of charged polarons. The formation of polaronic 
^ ' droplets in this case is due to competition of two types of interactions: the 

long range Coulomb repulsion and attraction generated by the deformation 
field. It is important to underline that if we consider the system of neutral po- 
larons (without Coulomb repulsion), it shows a first order phase transition at 
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constant chemical potential, and is unstable with respect to global phase sep- 
aration at fixed density [T7j. Electron-phonon interaction may be short range 
and long range, depending on the type of phonons involved. In most of the 
cases we consider phonons of the molecular type, leading to short range forces. 
In some cases we consider long-range Frohlich electron-phonon interaction and 
interaction with the strain. 



2 Single polaron in the adiabatic approximation 



Adiabatic theory of polarons was formulated many years ago [IHl [19] and 
we briefly formulate here the main principles of adiabatic theory for the par- 
ticular case of interaction with molecular vibrations (Holstein model). The 
central equation in the adiabatic theory is the Schrodinger equation for the 
electron in the external potential of the deformation field. In the discrete 
version it has the following form pUj : 

- J2 ti^n^'r. - ^i+J + V2gu^^J,i = Eu^l (1) 

Here <(m) is hopping integral, V'n is electronic wave function on the site n, 
</5n is the deformation at the site n, g is electron phonon coupling constant 
and Wq is phonon frequency, k describes quantum numbers of the problem. 
Important assumption of the adiabatic approximation is that deformation 
field is very slow and we assume that (/?n is time independent d^p/dt = when 
we substitute it to the Schrodinger equation for electron. Therefore, — > 
and g'^ oo but the product g'^ujQ — Ep is finite and called polaron shift. The 
equation for (fn has the form [20j : 

= (2) 

Here corresponds to the ground state solution of the Eq.(l). After substi- 
tution of the Eq.(2) to Eq.(l) we obtain: 

- ^ t(m)[V^^ - - 2£;,|^°|V^ = Eki^l (3) 

As a result nonlinear Schrodinger equation (Eq.(3)) describes the ground state 
of the polaron. All excited states are the eigenvalues and eigenfunctions of the 
linear Schrodinger equation in the presence of the external field determined 
by the deformation field Eq.(2). Polaron energy is the sum of two contribu- 
tions. The first is the energy of electron in the self-consistent potential well, 
determined by Eq.(3), and the second is the energy of the strain field ip itself 
Epoi = Eq + The polaron energy is presented in the Fig.l as a 



Polarons; from single polaron to short scale phase separation. 



3 




Fig. 1. Polaron energy as a function of coupling constant in ID, 2D and 3D cases. 
Dashed lines represent the energy of delocalized solution in ID, 2D and 3D respec- 
tively 



function of the dimcnsionlcss coupling constant Ig^uji^lt in ID, 2D and 3D 
cases. 

There is very important difference between ID and 2, 3D cases. Polaron en- 
ergy for ID case is always less then the energy of the delocalized state(dashed 
line). Polaron is always stable in ID. In 2D and 3D cases there is critical value 
of the coupling constant where first localized solution of Eq.(3) appears. It is 
interesting that the energy of the solution is higher then energy of the delocal- 
ized state. Therefore in 2D and 3D cases there is range of coupling constant 
where polaron is metastable. Delocalized solution is always stable in 3D case. 
Therefore, the barrier, which separates localized and delocalized states exists 
in the whole region of the coupling constants, where self-trapped solution ex- 
ists. In 2D case the delocalized state is unstable at large value of the coupling 
constant. The barrier separating localized and delocalized states forms only 
in the restricted region of the coupling constant [20] . To demonstrate that we 
have plotted polaron energy as a function of its radius in 2D in the Fig. 2. As it 
is clearly seen from this figure, barrier has disappeared at g > gcz — \/2'Kt/ijjQ 

m- 

There are two types of non-adiabatic corrections to adiabatic polaron. 
The first is related to renormalization of local phonon modes. Fast motion 
of the electron within polaronic potential well leads to the shift of the local 
vibrational frequency: 

u: = uoa[l - ze /2{g^uo^ff'^ (4) 

here z is the number of nearest neighbors. This formula is valid in the strong 
coupling limit g^ujQ ^ t and when tunneling frequency of the polaron is much 
smaller then phonon frequency. 
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Fig. 2. Polaron energy as a function of radius in 2D. For g > dEpoi/dR < 0. 
gel = 1.69i/t/a;o, and gc2 = 1.87i/I7^ 



Another type of correction is related to the restoration of translational 
symmetry (Goldstone mode) and describes polaron tunnelling and forma- 
tion of the polaron band. This correction was calculated in the original pa- 
per of Holstein (Ref . ^18j ) . Slightly improved formula was derived in |21j . In 
the adiabatic limit polaron tunnelling is exponentially suppressed teff oc 
EpUjQ exp {—g^) (see Eq. (9) of the Ref. [21]) and should be smaller then 
phonon frequency wo- 

In the following we will neglect all nonadiabatic corrections. We will con- 
sider polaron as pure localized state, and all corrections which contain phonon 
frequency itself and tunnelling amplitude for polaron are neglected. 



3 Strings in charge-transfer Mott insulators: effects of 
lattice vibrations and the Coulomb interaction 

Here we prove that the Frohlich electron-phonon interaction combined with 
the direct Coulomb repulsion does not lead to charge segregation like strings 
in doped narrow-band insulators, both in the nonadiabatic and adiabatic 
regimes. However, this interaction significantly reduces the Coulomb repul- 
sion, which might allow much weaker antiferromagnetic and/or short-range 
electron-phonon interactions to segregate charges in the doped insulators, as 
suggested by previous studies [H [51 [TT] . 

To begin with, we consider a generic Hamiltonian, including, respec- 
tively, the kinetic energy of carriers, the Frohlich electron-phonon interaction, 
phonon energy, and the Coulomb repulsion as 
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H = ^ t{m - n)6s^s'clcj + ^ uj^n, [ui(q)dq + H.c] 

+ '^q(4'^q + 1/2) + 2 II ^(m - (5) 

q i=ij 

with bare hopping integral t(m), and matrix element of the electron- phonon 
interaction 

u,iq) = ^^(q)e*q-. (6) 

Here i — (m, s), j — (n, s') include site m, n and spin s, s' quantum numbers, 
Hi = c|ci, Ci,dq are the electron (hole) and phonon operators, respectively, 
and N is the number of sites. At large distances ( or small q) one finds 

7(q)^q = — ^, (7) 

and 

V{m~n) = —^ 1- (8) 

eoo|m- n| 

The phonon frequency lu^ and the static and high-frequency dielectric con- 
stants in = e'^ — ^ are those of the host insulator {H = c = 1). 

In the adiabatic limit one can apply a discrete version of the continuos 
nonlinear equation [52] proposed for the Holstein model (Eq.(l)), extended 
to the case of the deformation and Frohlich interactions in Ref. [TTl ^ [TU] . 
Applying the Hartree approximation for the Coulomb repulsion, the single- 
particle wave-function, -!/;„ (the amplitude of the Wannier state |n)) obeys the 
following equation 

- Y t{m)[lpn ~ Ipn+m] - e(j)n1pn = Elpn- (9) 
m=£0 

The potential (j>n,k acting on a fermion k at the site n is created by the 
polarization of the lattice f. and by the Coulomb repulsion with the other 
M — 1 fermions, (j)'^ ^, 

Both potentials satisfy the descrete Poisson equation as 

M 

KZ\4fc=47re^|V'„,p|', (11) 
p=i 

and 

M 

eoo^0^,fc - -4^e J2 \^n,p\^ (12) 
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with Zi0„ — X]m('?^n ^ 0n+m)- Differently from Ref. ^ we include the 
Coulomb interaction in Pekar's functional J [55], describing the total energy, 
in a selfconsistent manner using the Hartree approximation, so that [10] 



J=- ^ Cp*(m)[-i/'n,p - V'n+m,p] 



27re2 



n,p,m,g 

+ ^ E IV'n.p|M-l|V'™,,p. (13) 

If we assume, following Ref. [TT] that the single-particle function of a 
fermion trapped in a string of the length is a simple exponent, V'n = 
A'^""'^/^ exp(ifcn) with the periodic boundary conditions, then the functional J 
is expressed as J = T+C/, where T = -2i(iV- 1) sin(7rM/7V)/[Af sin(7r/7V)] is 
the kinetic energy (for an odd number M of spinless fermions) , proportional 
to and 

[/ ^ Af2/^ + — M{M - 1)In, (14) 

corresponds to the polarisation and the Coulomb energies. Here the integral 
In is given by 

— . o / dx [ dy I dz ' ^ ^ ^ 



(2^)3 _iV2sin(:r/2)2 
X (3 — cosx — cosy — cosz) (15) 

Im has the following asymptotic [T^: 

1.31 + In iV 

-^w = > (16) 

The asymptotic is derived also analyticially at large N by the use of the 
fact that sin(7Vx/2)2/(27riVsin(x/2)2) can be replaced by a 8- function. If we 
split the first (attractive) term in Eq.(14) into two parts by replacing AP for 
M + M{M — 1), it becomes clear that the net interaction between polarons 
remains repulsive in the adiabatic regime because k > eoo- Hence, there are no 
strings within the Hartree approximation for the Coulomb interaction. Strong 
correlations do not change this conclusion. Indeed, if we take the Coulomb 
energy of spinless one-dimensional fermions comprising both Hartree and ex- 
change terms as 

Ec = ^''^y^ ~ [0.916 + InM], (17) 
the polarisation and Coulomb energy per particle becomes (for large M >> 1) 

U/M = -^-^[0.916 + InM - a(l. 31 + ln7V)], (18) 
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where a — I — £00/^0 < 1- Minimising this energy with respect to the length 
of the string N we find 

iV = Mi/"exp(-0.3f + 0.9f6/a), (f9) 

and 

g2 

{U/M)rmn = M^-^/" exp(0.3f - 0.9f 6/a). (20) 

K 

Hence, the potential energy per particle increases with the number of particles 
so that the energy of M well separated polarons is lower than the energy of 
polarons trapped in a string no matter whether they are correlated or not. 
The opposite conclusion of Ref. [9j originates in an incorrect approximation of 
the integral In oc N'-'-^^/N. The correct asymptotic result is In — \n{N)/N. 

One can argue [9 that a finite kinetic energy (t) can stabilise a string of 
a finite length. Unfortunately, this is not correct either. We performed exact 
(numerical) calculations of the total energy E{M, N) of M spinless fermions in 
a string of the length N including both kinetic and potential energy with the 
typical values of e^o = 5 and eo = 30. The local energy minima (per particle) 
in the string of the length 1 < < 69 containing M < N/2 particles are 
presented in the Table 1. Strings with even fermion numbers carry a finite 
current and hence the local minima are found for odd M. In the extreme 
wide band regime with t as large as 1 eV the global string energy minimum 
is found at M = 3,iV = 25 {E = -2.1167 eV), and at Af = 3,iV = 13 for 
t = 0.5 eV {E = —1.2138 eV). However, this is not the ground state energy in 
both cases. The energy of well separated d > 2-dimensional polarons is well 
below, less than —2dt per particle (i.e. —6 eV in the first case and —3 eV in 
the second one in the three dimensional cubic lattice, and —4 eV and —2 eV, 
respectively, in the two-dimensional square lattice). This argument is applied 
for any values of eoi^oo and t. As a result we have proven that strings are 
impossible with the Frohlich interaction alone contrary to the erroneous Ref. 



Table 1. E{M, N) for t = 1 eV and t = 0.5 eV 





t = 


leV 


t = 


0.5eV 


M 


N 


E(M,N) 


N 


E(M,N) 


1 


11 


-2.0328 


3 


-1.1919 


3 


25 


-2.1167 


13 


-1.2138 


5 


42 


-2.1166 


25 


-1.1840 


7 


61 


-2.1127 


40 


-1.1661 



The Frohlich interaction is, of course, not the only electron-phonon interac- 
tion in ionic solids. As discussed in Ref. p3], any short range electron-phonon 
interaction, like, for example, the Jahn- Teller (JT) distortion can overcome 
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the residual weak repulsion of Frohlich polarons to form small bipolarons. 
At large distances small nonadibatic bipolarons weakly repel each other due 
to the long-range Coulomb interaction, four times of that of polarons, Eq.(9). 
Hence, they form a liquid state !^ , or bipolaronic-polaronic crystal-like struc- 
tures [24] depending on their effective mass and density. The fact, that the 
Frohlich interaction almost nullifies the Coulomb repulsion in oxides justifies 
the use of the Holstein-Hubbard model [25l [26] . The ground state of the ID 
Holstein-Hubbard model is a liquid of intersite bipolarons with a significantly 
reduced mass (compared with the on-site bipolaron) as shown recently [57]. 
The bound states of three or more polarons are not stable in this model, thus 
ruling out phase separation. However, the situation might be different if the 
antiferromagnetic [4[ [5] and JT interaction f6| or any other short (but finite) 
range electron-phonon interaction are strong enough. Due to long-range na- 
ture of the Coulomb repulsion the length of a string should be finite (see, also 
Ref.[II[[Tni[5S]). 

To summarize we conclude that there are no strings in ionic doped in- 
sulators with the Frohlich interaction alone. Depending on their density and 
mass polarons remain in a liquid state or Wigner crystal. On the other hand 
the short-range electron-phonon and/or antiferromagnetic interactions might 
provide a liquid bipolaronic state and/or charge segregation (strings of a fi- 
nite length) since the long-range Frohlich interaction significantly reduces the 
Coulomb repulsion in highly polarizable ionic insulators. 

4 Ordering of charged polarons: Lattice gas model 

In this section we consider macroscopic system of polarons in the thermo- 
dynamic limit. To underline nontrivial geometry of the phase separation we 
consider two-fold degenerate electronic states which interact with nonsym- 
metric deformation field. In our derivation we follow particular model for 
high-Tc superconductors. Nevertheless the results are general enough and are 
applicable to many Jahn- Teller systems. 

Recently we formulated the model[29] where we suggested that interac- 
tion of a two-fold degenerate electronic state with fully symmetric of the 
small group ri phonon modes at a finite wave-vector can lead to a local non- 
symmetric deformation and short-length scale charge segregation in high-Tc 
materials. Here we reduced the proposed model to the lattice gas model [T7] 
and showed that the model indeed displays phase separation, which may oc- 
cur in the form of stripes or clusters depending on the anisotropy of the 
short range attraction between localized carriers [l^. We also generalized the 
model taking into account interaction of the Jahn- Teller centers via elastic- 
ity induced field [17] . We showed that the model without Coulomb repulsion 
displays a first order phase transition at a constant chemical potential. When 
the number of particles is fixed, the system is unstable with respect to the 
global phase separation below certain critical temperature. In the presence of 
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the Coulomb repulsion global phase separation becomes unfavorable due to 
a large contribution to the energy from long range Coulomb interaction. The 
system shows mesoscopic phase separation where the size of charged regions 
is determined by the competition between the energy gain due to ordering 
and energy cost due to breaking of the local charge neutrality. Since the short 
range attraction is anisotropic the phase separation may be in the form of 
short segments or/and stripes. 

Let us start with the construction of a real-space Hamiltonian which cou- 
ples 2-fold degenerate electronic states (or near-degenerate states) with op- 
tical phonons of ti symmetry. Two-fold degeneracy is essential because in 
that case formation of the polaronic complexes leads to reduction not only 
of translational symmetry but also reduction of the point group symmetry. 
Since the Hamiltonian needs to describe a 2-fold degenerate system, the 2- 
fold degenerate states - for example the two £"„ states corresponding to the 
planar hybridized Cu d^^-y^ , O Px and Py orbitals, or the and Eg states 
of the apical O - are written in the form of Pauli matrices ai. Taking into 
account that the states are real, the Pauli matrices which describe transitions 
between the levels transform as Aig {k^ + ky) for cto, Big (fc^ — ky) for 0-3, 
B2g (kxky) for (Ti, and (sz) for (T2 representations respectively. Collecting 
terms together by symmetry we can construct effective electron-spin-lattice 
interaction Hamiltonian given in the Ref. [29]. Here we consider a simplified 
version of the JT model Hamiltonian [29] , taking only the deformation of the 
Big symmetry: 

HjT=9Y,a:,,,f{n){b[^^ + b,+^), (21) 

n,l 

here the Pauli matrix 173^1 describes two components of the electronic doublet, 
and /(n) = (n^ — n'^)fo(n) where f(i(n) is a symmetric function describing 
the range of the interaction. For simplicity we omit the spin index in the 
sum. The model could be easily reduced to a lattice gas model[17]. Let us 
introduce the classical variable 'Pi =< bf + hi> and minimize the energy 
as a function of in the presence of the harmonic term w ^? /2. We obtain 
the deformation, corresponding to the minimum of energy, 

'Pf^ = -V25/c^^a3,i+„/(n). (22) 

n 

Substituting <pf'^ into the Hamiltonian (1) and taking into account that the 
carriers are charged we arrive at the lattice gas model. We use a pseudospin 
operator = 1 to describe the occupancies of the two electronic levels niand 
712- Here S*^ = 1 corresponds to the state with ni = 1 , 77.2 = 0, S'f = — 1 
to ni = 0, 712 = 1 and 5f = to = 712 = 0. Simultaneous occupancy of 
both levels is excluded due to a high on-site Coulomb repulsion energy. The 
Hamiltonian in terms of the pseudospin operator is given bv[17) 

Hlg = Y.{-Vi{\-i)SlSI -f V,{i - miQi), (23) 
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where Qi = (5'f K;(n) = e^/eoa(n-^ + "y)^^^ is the Coulomb potential, e is 
the charge of electron, cq is the static dielectric constant and a is the effective 
unit cell period. The anisotropic short range attraction potential is given by 
V;(n) = ff^/t^ X^m /('^)/('^+™)- The attraction in this model is generated by 
the interaction of electrons with optical phonons. The radius of the attraction 
force is determined by the radius of the electron-phonon interaction and the 
dispersion of optical phonons [TP]. 

A similar model can be formulated in the limit of continuous media [T7]. 
The deformation is characterized by the components of the strain tensor. For 
the two dimensional case we can define 3 components of the strain tensor: 
ei — Uxx + Uyy, e — Uxx — Uyy and 62 — u^y transforming as the Aig, Big and 
-B25 representations of the Z?4/i group respectively. These components of the 
tensor are coupled linearly with the two-fold degenerate electronic state which 
transforms as the Eg or Eu representation of the point group. Similarly to the 
case of previously considered interaction with optical phonons we keep the 
interaction with the deformation e of the Big symmetry only. The Hamiltonian 
without the Coulomb repulsion term has the form: 

H^gY, Sfe, + i {Aiel, + A^e? + Age^^;) , (24) 

i 

where Aj are corresponding components of the elastic modulus tensor. The 
components of the strain tensor are not independent [SUl 113 and satisfy the 
compatibility condition: 

V2ei(r) - 4d^e2ir)/dxdy = {dydx^ ~ d"^ / dy'^)e{v) 

The compatibility condition leads to a long range anisotropic interaction be- 
tween polarons. The Hamiltonian in the reciprocal space has the form: 

el 



H^gY,Slei, + {A2+AiUm^. (25) 



The Fourrier transform of the potential is given by: 



By minimizing the energy with respect to and taking into account the long- 
range Coulomb repulsion we again derive Eq.(23). The anisotropic interaction 

2 

potential Vi (n) = — exp (ik • n) 2{A-2+AiU('k)) determined by the inter- 
action with the classical deformation and is long-range. It decays as at 
large distances in 2D. Since at large distances the attraction forces decay faster 
then the Coulomb repulsion forces the attraction can overcome the Coulomb 
repulsion at short distances, leading to a mesoscopic phase separation. 

Irrespective of whether the resulting interaction between polarons is gen- 
erated by acoustic or optical phonons the main physical picture remains the 
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same. In both cases there is an anisotropic attraction between polarons on 
short distances. The interaction could be either ferromagnetic or antiferromag- 
netic in terms of the pseudospin operators depending on mutual orientation 
of the orbitals. Without loosing generality we assume that ^(n) is nonzero 
only for the nearest neighbors. 

Our aim is to study the model (Eq.(23)) at constant average density, 



^E^i' (27) 



N 
1 

where N is the total number of sites. However, to clarify the physical picture 
it is more appropriate to perform calculations with a fixed chemical potential 
first by adding the term — Qi ^o the Hamiltonian (23). 

Models such as Eq.(23) without the long-range forces, were studied many 
years ago on the basis of the molecular-field approximation in the Bragg- 
Williams formalism (STJ [32] . The mean-field equations for the particle density 
n and the pseudospin magnetization M = J^i have the form[3Tj: 

M = 2sinh(2zV^A//fc,T) 

exp{~fi/kgT) + 2cosh{2zViM/kBT) ^ ' 



2cosh{2zViM/kBT) 



e-xp {-fi/ksT) + 2 cosh {2zViM/kBT) 



(29) 



here z = 4 is the number of the nearest neighbours for a square lattice in 
2D and ks is the Boltzman constant. A phase transition to an ordered state 
with finite M may be of either first or second order, depending on the value 
of fi. For the physically important case —2zVi < /i < 0, ordering occurs as 
a result of the first order phase transition. The two solutions of Eqs. (28,29) 
with M = and M 7^ correspond to two different minima of the free energy. 
The temperature of the phase transition Tcrit is determined by the condition: 
F{M = 0, /X, T) = F(M, T) where M is the solution of Eq. (28). When the 
number of particles is fixed (Eq.(29)), the system is unstable with respect to 
global phase separation below Tcrit- As a result, at fixed n two phases coexist 
with rto — n{M = 0, /i, T) and um = n(M, fi, T), resulting in a liquid-gas- like 
phase diagram (Fig. 3). 

To investigate the effects of the long range-forces, we performed Monte 
Carlo simulations on the system (23) [1^. The simulations were performed 
on a square lattice with dimensions up to L x L sites with 10 < L < 100 
using a standard Metropolis algorithm[33j in combination with simulated 
annealing [34]. At constant n one Monte Carlo step included a single update 
for each site with nonzero Qi, where the trial move consisted from setting 
5*2 = at the site with nonzero and Sz — ±lat a randomly selected site 
with zero Qi. A typical simulated annealing run consisted from a sequence of 
Monte Carlo simulations at different temperatures. At each temperature the 
equilibration phase (10"^ — 10^ Monte Carlo steps) was followed by the averag- 
ing phase with the same or greater number of Monte Carlo steps. Observables 
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were measured after each Monte Carlo step during the averaging phase only. 
For L > 20 we observe virtually no dependence of the results on the system 
size. 

Comparing the Monte Carlo results in the absence of Coulomb repulsion 

shown by tc-it in Fig- 3 with MF theory we find the usual reduction of tcrit 
due to fluctuations in 2D by a factor of ~ 2. 



I h 




0.0 0.2 0.4 0.6 0.8 1.0 1 2 
n t 

Fig. 3. a) The phase diagram generated by Hjt (23) with, and without the Coulomb 
repulsion (CR). The dashed line is the MF critical temperature, while the full tri- 
angles (a) represent the Monte Carlo critical temperature, tcru, without CR. The 
open circles (o) represent td, without CR. The open triangles (A) represent td while 
the diagonal crosses (x) represent the onset of clustering, to, in the presence of CR. 
The cluster-ordering temperature (see text), tco, (also incl. CR) is shown as crosses 
{+). The size of the symbols corresponds to the error bars, b) Typical temperature 
dependencies of the nearest neighbor density correlation function for n = 0.18 
in the absence of CR (•) and in the presence of CR (V). Arrows indicate the char- 
acteristic temperatures. 



Next, we include the Coulomb interaction Vc{r). We use open bound- 
ary cx)ii(liti()iis to avoid complications due to the long range Coulomb forces 
and ensure overall electroneutrality by adding a uniformly charged back- 
ground electrostatic potential (jellium) to Eq. (23). The short range potential 
Vi{i) ~ Vi{\)eQale? was taken to be nonzero only for |i| < 2 and is therefore 
specified only for nearest, and next-nearest neighbours as i'((l,0) and 1) 
respectively. 

The anisotropy of the short range potential has a profound influence on the 
particle ordering. We can see this if we fix vi{\, 0) = —1, at a density n = 0.2 
and vary the next-nearest neighbour potential fi(l, 1) in the range from — 1 
to 1. When W((l,l) < 0, the attraction is " ferrodistortive" in all directions, 
while for positive v;(l, 1) > the interaction is "antiferrodistortive" along the 
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diagonals. The resulting clustering and ordering of clusters att — 0.04 is shown 
in Fig. 4a). As expected, a more symmetric attraction potential leads to the 
formation of more symmetric clusters. On the other hand, for vi{l, 1) = 1, the 
" antiferrodistortive" interaction along diagonals prevails, resulting in diagonal 
stripes. 

In the temperature region where clusters partially order the heat capac- 
ity (cL = d{E)L/dT where E is the total energy) displays the peak at tco- 
The peak displays no scaling with L indicating that no long range ordering 
of clusters appears. Inspection of the particle distribution snapshots at low 
temperatures (Fig. 4a) reveals that finite size domains form. Within domains 
the clusters are perfectly ordered. The domain wall dynamics seems to be 
much slower than our Monte Carlo simulation timescale preventing domains 
to grow. The effective L is therefore limited by the domain size. This explains 
the absence of the scaling and clear evidence for a phase transition near tco- 

We now focus on the shape of the short range potential which pro- 
motes the formation of stripes shown in Fig. 4a). We set vi{l,0) = — 1 
and vi{l,l) = and study the density dependence. Since the inclusion of 
the Coulomb interaction completely suppresses the first order phase tran- 
sition at tcriti we measure the nearest neighbor density correlation function 
9pL = 4n(i-n)L^ E|m|=i (Ei (Qi+m - n) {Qi - n)) ^ to dctcct clustcriug. Here 
0^ represents the Monte Carlo average. We define a dimensionless temper- 
ature tci = ksTcieoa/e^ as the characteristic crossover temperature related 
to the formation of clusters at which gpL rises to 50% of its low temperature 
value. The dependence of td on the density n is shown in the phase diagram 
in Fig. 3. Without Coulomb repulsion Vc{r), td follows tcrit, as expected. The 
addition of Coulomb repulsion Vc{r) results in a significant decrease of td 
and suppression of clustering. At low densities we can estimate the onset for 
cluster formation by the temperature, ip: at which gpj^ becomes positive. It is 
interesting to note that t^ almost coincides with the tcrit line at low n (Fig. 
3). 

To illustrate this behaviour, in Fig. 4b) we show snapshots of the cal- 
culated Monte Carlo particle distributions at two different temperatures for 
different densities. The growth and ordering of clusters with decreasing tem- 
perature is clearly observed. At low n, the particles form mostly pairs with 
some short stripes. With further increasing density, quadruples gradually re- 
place pairs, then longer stripes appear, mixed with quadruples, etc.. At the 
highest density, stripes prevail forming a labyrinth-like pattern. The density 
correlation function shows that the correlation length increases with doping, 
but long range order is never achieved (in contrast to the case without Vc). 
Note that while locally there is no four-fold symmetry the overall correlation 
function still retains 4-fold symmetry. 

To get further insight in the cluster formation we measured the cluster-size 
distribution [IT]. In Fig. 5 we present the temperature and density dependence 
of the cluster-size distribution function xlU) = {Np{j)) ^ / [nL"^), where Np{i) 
is the total number of particles within clusters of size j. At the highest tem- 
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a) 



b) /) = 0.08 n = 0.26 

t = 064 .■ -^i^r^^'^^r " 



Fig. 4. a) Snapshots of clusters ordering at t — 0.04, n — 0.2 and w;(l,0) = — 1 
for different diagonal vi{l, 1) (given in each figure). Grey and black dots represent 
particles clusters in state Si = 1 and states Si = — 1 respectively. The preference 
for even-particle-number clusters in certain cases is clearly observed, for example 
for vi{l, 1) = —0.2. b) Snapshots of the particle distribution for two densities at two 
different temperatures t — 0.64 and t — 0.1 respectively. 

perature XlU) is close to the distribution expected for the random ordering. 
As the temperature is decreased, the number of larger clusters starts to in- 
crease at the expense of single particles. Remarkably, as the temperature is 
further reduced, clusters of certain size start to prevail. This is clearly seen 
at higher densities (Fig. 3). Depending on the density, the prevailing clusters 
are pairs up to n w 0.2, quadruples for 0.1 ^ n < 0.3 etc.. We note that 
for a large range of w;(l,0), the system prefers clusters with an even number 
of particles. Odd particle-number clusters can also form, but have a much 
narrower parameter range of stability. The preference to certain cluster sizes 
becomes clearly apparent only at temperatures lower then td, and the tran- 
sition is not abrupt but gradual with the decreasing temperature. Similarly, 
with increasing density changes in textures also indicate a series of crossovers. 

The results of the Monte Carlo simulation[T7] presented above allow a 
quite general interpretation in terms of the kinetics of first order phase 
transitions [35]. Let us assume that a single cluster of ordered phase with 
radius R appears. As was discussed in [71 [29], the energy of the cluster is 
determined by three terms: e = —FnE? + anR + ^R^ . The first term is the 
energy gain due to the ordering phase transition where F is the energy dif- 
ference between the two minima in the free energy density. The second term 
is the surface energy parameterized by a, and the third term is the Coulomb 
energy, parameterized by 7. If a < 7tF/3j, e has a well defined minimum at 
R = Rq corresponding to the optimal size of clusters in the system. Of course, 
these clusters are also interacting among themselves via Coulomb and strain 
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f n 

Fig. 5. The temperature dependence of the cluster-size distribution function XL{j) 
(for the smallest cluster sizes) as a function of temperature at two different average 
densities n — 0.08 (a) and n = 0.18 (b). XL{j) as a function of n at the temperature 
between to and td (c), and near tco (d). The ranges of the density where pairs prevail 
are very clearly seen in (d). Error bars represent the standard deviation. 

forces, which leads to cluster ordering or freezing of cluster motion at low 
temperatures as shown by the Monte Carlo simulations. 

We conclude that a model with only anisotropic JT strain and a long- 
range Coulomb interaction indeed is unstable with respect to the short scale 
phase separation and gives rise to a remarkably rich phase diagram includ- 
ing pairs, stripes and charge- and orbital- ordered phases, of clear relevance 
to oxides. The energy scale of the phenomena is defined by the parameters 
used in Hjt-c (23). For example, using the measured value eg — 40 |35] for 
La2Cu04, we estimate 0) = 0.1 eV, which is also the typical energy scale 
of the "pseudogap" in the cuprates. The robust prevalence of the paired state 
in a wide region of parameters (Fig. 5 c,d) is particularly interesting from 
the point of view of superconductivity. A similar situation occurs in mangan- 
ites and other oxides with the onset of a conductive state at the threshold 
of percolation, but different textures are expected to arise due to different 
magnitude (and anisotropy) of V;(n), and static dielectric constant cq in the 
different materials |37j. 

5 Coulomb frustrated first order phase transition 

As it is stated in the previous section uncharged JT polarons have the ten- 
dency to ordering. The ordering transition is a phase transition of the first or- 
der. At the fixed density of polarons the system is unstable with respect to the 
global phase separation. The global phase separation is frustrated by charg- 
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ing effects leading to short-scale phase separation. Therefore the results of the 

Monte-Carlo simulation of the model (23) allow general model independent 
interpretation. Let us consider the classical free energy density corresponding 
to the first order phase transition: 

F, = {{t - 1) + (ry2 _ 1)2)^2 (30) 

Here t = {T - Tc)/{Tq - T^) is dimensionless temperature. At t = 4/3 (T = 
To + (To — T'c)/3) the nontrivial minimum in the free energy appears. At 
t = l{T = To) the first order phase transition occurs, but the trivial solution 
77 = corresponds to the metastable phase. At t = (T = Tc) trivial solution 
becomes unstable. In order to study the case of Coulomb frustrated phase 
transition we have to add coupling of the order parameter to local charge 
density. Our order parameter describes sublattice magnetization and therefore 
only square of the order parameter may be coupled to the local charge density 

Fcoupi = -aifp (31) 

The total free energy density should contain the gradient term and the elec- 
trostatic energy: 

Fgrad + Fei = C{Vvf + f [p(r) - p] J dr' [p{r' ) - p]/\r - r'| (32) 

Here we write p explicitly to take into account global electroneutrality. Total 
free energy Eqs. (30-32) should be minimized at fixed t and p. 

Let us demonstrate that Coulomb term leads to the phase separation in 
the 2D case. Minimization of F with respect to the charge density p(r) leads 
to the following equation: 

- aVlnV^ = 4nK[p{r) - p]d5{z) (33) 

here we write explicitly that density p(r) depends only on 2D vector r and 
introduce layer thickness d, to preserve correct dimensionality. Solving this 
equation by applying the Furrier transform and substituting the solution back 
to the free energy density we obtain: 

. = ..-o,V.W-^/..«p (34, 

As a results the free energy functional is similar to the case of the first order 
phase transition with shifted critical temperature due to the presence of the 
term arj'^p and nonlocal gradient term of higher order. 

To demonstrate that uniform solution has higher energy then nonhomo- 
geneous solution we make Fourier transform of the gradient term: 

F,™dOcCfc^hkp-^^^^j^^ (35) 
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here r/k a-nd (?7^)k are are Fourier components of the order parameter and 
square of the order parameter respectively. If we assume that the solution is 
uniform i.e. Typ 7^ and (??^)o 7^ small nonuniform corrections to the solution 
reduces the free energy at small k, where second term dominates. 

Proposed free energy functional is similar to that proposed in the Ref. [16| . 
The important difference is that in our case the charge is coupled to the square 
of the order parameter and it plays the role of local temperature, while in the 
case of Ref.[16] there is linear coupling of the charge to the order parameter. 
Charge in that case plays the role of the external field. Moreover, contrary to 
the case of Ref[T6] where charge is accumulated near domain walls, in our 
case charge is accumulated near interphase boundaries. 

6 Conclusion 

We have demonstrated that anisotropic interaction between Jahn- Teller cen- 
ters generated by optical and/or acoustical phonons leads in the presence of 
the long range Coulomb repulsion to the short scale phase separation. Topol- 
ogy of texturing differs from charged bubbles to oriented charged stripes de- 
pending on the anisotropy of short range potential. On the phenomenological 
level inhomogeneous phase with charged regions appears due to tendency of 
the system of polarons to global phase separation, while global phase sepa- 
ration is frustrated by the long-range Coulomb forces. Effectively this system 
may be described by standard Landau functional with nonlocal long-range 
gradient term. 
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